Quantum Dynamics of Solitons in Strongly Interacting Systems on Optical Lattices 
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Mean-field dynamics of strongly interacting bosons has been shown to support two species of solitons: one of 
Gross-Pitaevski (GP)-type where the condensate fraction remains dark and a novel non-GP-type characterized 
by brightening of the condensate fraction. Here we study the effects of quantum fluctuations on these solitons 
using the adaptive time-dependent density matrix renormalization group method, which takes into account the 
effect of strong correlations. We use local observables as the density, condensate density and correlation func- 
tions as well as the entanglement entropy to characterize the stability of the initial states. We find both species of 
solitons to be stable under quantum evolution for a finite duration, their tolerance to quantum fluctuations being 
enhanced as the width of the soliton increases. We describe possible experimental realizations in atomic Bose 
Einstein Condensates, polarized degenerate Fermi gases, and in systems of polar molecules on optical lattices. 



PACS numbers: 03.75.Lm, 03.75.-b, 67.85.De 

I. INTRODUCTION 

Solitary waves and solitons (i.e., solitary waves whose 
shape remains unchanged at collisions) are encountered in 
systems as diverse as classical water waves [1], magnetic 
materials [1-17], fiber-optic communication [1, 18, 19], and 
Bose-Einstein condensates (BEG) [1, 20-22]. Rooted in the 
nonlinearity of the system which balances dispersive effects, 
solitons are fascinating non-linear waves that encode collec- 
tive behavior in the system. Intrinsically nonlinear in nature 
due to inter-particle interactions, and due to the high degree 
of control in the experiments, the BEG systems are a natural 
fertile ground for exploring solitons. In dilute atomic gaseous 
BECs which are simply described in terms of the properties 
of the non-linear Schrodinger equation (NLSE), or Gross- 
Pitaevski equation (GPE), bright (density elevation) solitons 
exist for attractive interparticle interactions and dark (density 
notch) solitons in the repulsive case. These solutions are char- 
acterized not only by persistent density profiles, but also by 
characteristic modulations of the quantum phase across their 
profiles, which differs for the bright (attractive condensate) 
[23] and the dark (repulsive condensate) cases [24-26]. How- 
ever, we want to emphasize that, on general grounds, other 
systems with intrinsic nonlinearities should be able to realize 
solitons if the conditions are chosen appropriately. 

In this paper, we follow two goals: First, we want to de- 
scribe the realization of solitons in lattice systems in a broader 
framework which includes BEG, quantum degenerate Fermi 
gases, hard-core bosons, and polar molecules on optical lat- 
tices, as well as certain condensed matter systems. The com- 
mon aspect of these various systems is that the soliton dy- 
namics can be described in terms of a simple S = 1/2 spin 
chain which, in turn, can be the effective model for a variety 
of situations, as the ones mentioned above. This description is 
footed on an extension of the standard GPE treatment of soli- 



tons, and leads us to the second scope of our paper which is to 
describe new effects which go beyond mean-field dynamics. 

Investigations of effects beyond GPE dynamics have been 
a subject of various studies in the past decade. For short range 
repulsive systems in the Tonks-Girardeau regime, the cubic 
non-linearity of GPE was replaced by a quintic repulsive non- 
linearity and the resulting modified GPE was shown to sup- 
port dark solitary waves of GP-type [27]. This ID system 
was further investigated in the presence of dipolar interactions 
[28] and was shown to support bright solitons whose stability 
and mobility depended on the dipolar interaction strength. In 
2D systems, bright solitons were found to be stable given a 
sufficient dipole-dipole strength [64]. Further studies of the 
stability and dynamics of solitons have been extended to two 
component BECs [65] and multilayered BECs [66]. Existence 
of dark and bright solitary waves was also shown numerically 
in systems describing multicomponent BECs [29]. In addi- 
tion to the continuum systems, solitary waves have been ex- 
tensively studied in systems described by a discrete non-linear 
Schrodinger equation (DNLSE) [30, 67], BECs in deep opti- 
cal lattices and also to optical beams in wave guides. 

In recent studies [47, 48], solitary waves in a system of hard 
core bosons (HCB) have been studied using mean field equa- 
tions obtained from mapping the HCB system to a 5* = 1/2 
spin chain and in 2D. The continuum limit of the lattice pop- 
ulated with HCB is described by a generalized GPE, which 
we will refer to as "HGPE" [See Eq. (6)]. In contrast to GPE, 
HGPE contains both the normal and condensate density. This 
system describing strongly repulsive BEC was shown to sup- 
port both dark and bright solitary waves, the existence of both 
species being rooted in the particle-hole duality/symmetry in 
HCB systems. Unlike other studies, HGPE solitary waves 
are obtained as an analytic solution which was shown to pro- 
vide an almost exact solution of the equations of motion [47] . 
These two species of solitons can be referred to as the GP- 
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type and the non-GP-type as the former corresponds to a dark 
condensate fraction that dies beyond sound velocity while the 
later is associated with brightening of the condensate and per- 
sists all the way up to sound velocity and transforms into a 
soliton train for supersonic velocities. Recent numerical stud- 
ies investigating collision properties of these nonlinear modes 
suggest that these solitary waves are in fact solitons [31]. 

An important question that we investigate here is whether 
these mean field solitons survive quantum fluctuations. In pre- 
vious work, the quantum dynamics of GP dark solitons in the 
superfluid regime of the Bose-Hubbard Hamiltonian (BHH) 
has been studied numerically by Mishmash et al. [32, 33]. 
The main findings are that for weak interactions the dark soli- 
ton is stable on a time-scale of the order of ~ 20 — 40 units of 
the hopping and afterwards decays due to two-particle scatter- 
ing processes. The larger the on-site interaction, the stronger 
the scattering and the faster the decay of the solitons. In ad- 
dition, these studies treated collisions between the solitary 
waves which confirm the soliton nature of the states on the 
time scales treated. These studies focus on the limit of small 
interactions. Here, we treat the strong coupling case and study 
the fate of the soliton solutions obtained in the HGPE frame- 
work. We do this by generalizing the BHH of Refs. [32, 33] 
to include on-site and nearest neighbor density-density inter- 
actions. As discussed in Ref. [47], in the continuum limit 
this gives rise to the two distinct types of solitons mentioned 
above, which, as we shall see, are found to be stable in both 
the mean field approximation to the lattice dynamics of the 
system, as well as in the full quantum dynamics on the lattice. 

More specifically, we describe the exact quantum evolu- 
tion of an initial mean field soliton solution on ID lattice sys- 
tems. The soliton and the Hamiltonian driving the dynamics 
are thereby formulated in terms of a 6* = 1/2 spin language. 
It is so possible to envisage a realization of the described soli- 
ton solutions in both, experiments with ultracold bosonic and 
spin polarized fermionic atoms, as well as in experiments with 
polar molecules [34-37] on optical lattices which can be used 
to emulate spin- 1/2 systems [38, 39]. We combine an analytic 
solution of the HGPE which provides a continuum approx- 
imation to the lattice problem, a numerical treatment of the 
mean-field equations on the lattice, and a full quantum treat- 
ment of the dynamics by applying the time-dependent DMRG 
[40-44]. Both, mean-field and numerical results indicate that 
for a certain range of parameters the solutions found are in- 
deed stable solitons on the time scale treated. The non-GP- 
type soliton is found to be somewhat less tolerant of quantum 
effects compared to the GP-type. We characterize the stabil- 
ity of the solitons by considering the entanglement in the sys- 
tem: since in our set-up the initial soliton solutions are product 
states on the lattice, the entanglement entropy [45] should re- 
main zero for a stable soliton solution and hence is a measure 
for the stability of the soliton in the course of the time evolu- 
tion. In addition, we consider correlation functions which, on 
similar grounds, can be used to characterize its stability. 

The paper is organized as follows. In section II, we intro- 
duce the effective spin model and its derivation from HCB and 
spinless fermions on a lattice, the dynamical equation (HGPE) 
that describes the continuum approximation to the mean-field 



equations of the lattice system, and we summarize the analytic 
solution of the HGPE. In Sec. Ill we describe the mean-field 
ansatz and some details of the DMRG approach to the dynam- 
ics. In Sec. IV, we analyze the stability of the soliton solutions 
by comparing the mean-field results on a lattice to the DMRG 
results. As a measure for the quality of the soliton solution, we 
use in Sec. IV B the von Neumann or entanglement entropy as 
well as correlation functions which also should remain zero in 
the course of the time evolution if the mean field state were to 
survive quantum fluctuations. In Sec. V we propose possible 
experimental realizations of the HGPE solitons. In Sec. VI, 
we summarize. 



II. HAMILTONIAN AND EQUATIONS OF MOTION 

In this paper, we treat the dynamics of initial soliton states 
driven by the spin Hamiltonian 



(1) 



on a one-dimensional lattice, i.e., we are treating the dynam- 
ics of a XXZ-chain with a global external magnetic field of 
magnitude g. One way to obtain this effective Hamiltonian is 
as the limiting case of the extended Bose Hubbard model. 
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Here, 6^^^ are the annihilation (creation) operators for a bo- 
son at the lattice site j, rij is the number operator, J/2 is the 
n.n. tunneling strength, and /i is the chemical potential. An 
attractive nearest-neighbor interaction V < is introduced to 
soften the effect of a strong onsite interaction \U\ 0. The 
HCB limit (\U\ ^ oo) corresponds to the constraint that two 
bosons cannot occupy the same site. This HCB system can 
then be mapped to the model Eq. (1) [51], where the two spin 
states correspond to two allowed boson number states |0) and 
|1), and setting g = J — V. Note that this is in contrast to 
the study of Refs. [32, 33] in which the quantum dynamics 
was investigated in the original Bose-Hubbard system and not 
for the effective model Eq. (1). This is interesting since the 
existence of the proposed soliton solutions for this spin model 
has implications for further systems than the ultracold bosonic 
atoms usually considered when describing soliton phenomena 
in cold gases. In particular it should be noted that the XXZ 
model in ID can be obtained using the Jordan- Wigner trans- 
form from a system of spinless fermions 
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with the fermionic annihilation (creation) operators on 
site j, and rij = c^-c- the density on site j. Therefore, it should 



(3) 



3 



be possible to investigate the soliton dynamics in experiments 
with bosonic atoms, in spin systems, and in fermionic sys- 
tems. In Sec. V we discuss possible implementations in ex- 
periments with cold gases. Note that both models, Eq. (1) 
and Eq. (3) are fundamental models for describing condensed 
matter systems such as quantum magnets and systems of itin- 
erant electrons. It is therefore conceivable that, in principle, 
the proposed soliton solutions can be realized in such systems 
as well. 

For simplicity, we set up our discussion in the framework 
of bosonic systems, without loosing generality. Then, the spin 
flip operators S"^ = Sx ^ iSy correspond to the annihilation 
and creation operators of the corresponding bosonic Hamil- 
tonian, bj Sj'. Thus the order parameter that describes 
a BEC wave function is tpj = {Sj~), where the expectation 
value is obtained using spin coherent states [51]. In this mean- 
field description, the evolution equation for the order parame- 
ter is obtained by taking the spin-coherent state average of the 
Heisenberg equation of motion. The spin coherent state \rj) 
at each site j can be parametrized as: 
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With this choice, the HCB system is mapped to 
a system of classical spins [47, 51] via S = 
(I sin((9) cos(^), I sm{0) sin(^), | cos(l9)). Note that 
the particle density pj and the condensate density pj satisfy 
the relation pj = pj p^, with p^ = I — pj the hole density. 
In this representation, tpj = ^/pje*"^. This mean field 
treatment is contrasted to the standard GPE derived from 
the Bose-Hubbard model by taking the expectation value of 
the Heisenberg equation of motion with Glauber coherent 
states [46]. We cast the equations of motion in terms of the 
canonical variables and Sj = cos{6j) = {1 — 2pj) and 
obtain 
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A. Solitary Waves in the Continuum Approximation 

In the continuum approximation, the equations for the or- 
der parameter are derived from a Taylor series in the lattice 
spacing a [54], 



ihi)' = (l-2pW^tP'-Vetp'V^p^UePtp'-ptp' (6) 
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where Ja^ = ^,Ue = 2(J - V) and K = Va'^. 

These equations have been shown to support solitary waves 







FIG. 1 . (Color online) Bright (top) and dark (bottom) soliton solu- 
tion in the continuum [Eqs. (7) and (9)]. Left panels show the den- 
sity, right panels the condensate density as a function of position and 
speed. Note that the condensate density of the bright soliton shows a 
'brightening', around the notch, i.e., it grows above the background 
value, whereas the dark soliton does not show this effect. 



[47] riding upon a background density po- p(^) = Po 
with z = X — vt. We obtain for the soliton solution 
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where 7 = Vl — a nd v being the speed of the solitary 
wave in units of Cg = y^2pq(1 — V/ J), which is the speed of 
sound of the Bose gas system determined from its Bogoliubov 
spectrum [51]. F is the width of the soliton. 
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The characteristic phase jump associated with the solitary 
waves is 



A(^± = (Vl - 2c2) cos- 
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1 - 2pg^2 • 
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This solution has some remarkable properties. One direct 
consequence of the particle-hole symmetry underlying the 
equations of motion is the presence of two species of solitary 
waves, shown in Fig. 1. The existence of /(z, po) superposed 
on the background particle density po implies the existence 
of a counterpart /(x,Pq), superposed upon a corresponding 
hole density p^. In fact it is easy to see that /^(z,po) = 
Po). For Po < 1/2, the ± corresponds to bright and 
dark solitons, respectively. The bright solitons have the un- 
usual property of persisting at speeds up to the speed of sound, 
in sharp contrast to the dark species that resembles the dark 
soliton of the GPE whose amplitude goes to zero at sound 
velocity. In the special case with background density equal to 
1/2, the two species of solitons become mirror images of each 
other, as/+(z,po = 1/2) = —f~{z^po = 1/2). In this case, 
the condensate density in fact describes the GPE-soliton [49]. 
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FIG. 2. (Color online) Comparison of the soliton profiles at times 
t — 20/ J obtained in the continuum [black solid line, Eq. (6)] and 
using the equations of motion approach Eq. (5) on a lattice of L = 
40 sites. The left panels show the results for V/ J — 0.75 (narrow 
soliton), the right panels the case V j J — 0.95 (broad soliton). The 
top panels show the bright soliton, the bottom panels the dark soliton 
solution of Eqs. (7) and (9). Both, the particle density p (*) and the 
condensate density />*(+) are shown. 



It should be noted that for po > 1/2, the dark and bright 
solitons switch their roles. In other words, for po < 1/2, it 
is the dark soliton that behaves like a GPE soliton while the 
bright soliton is the new type of soliton that persists all the 
way up to sound velocity. In contrast, for po > 1/2, the bright 
soliton is GPE-type while the dark one is the persistent soliton. 
In view of the particle-hole duality, we will present our results 
for po < 1/2 in which the bright solitons have the persistent 
character noted above. 

In the following sections, we will investigate for the ex- 
istence and the lifetime of these solutions on lattice systems 
using mean-field equations and the DMRG. We will comple- 
ment this analysis by investigating the stability of further ini- 
tial states. In particular, we show that an initial Gaussian den- 
sity distribution for a stationary soliton shows a similar stabil- 
ity if a phase jump is realized, but becomes unstable without 
a phase jump. This is of importance for experimental realiza- 
tions indicating that imperfections in the creation of the initial 
state may not have a strong influence on the soliton dynamics. 



III. METHODS: MEAN FIELD ANSATZ AND DMRG 



A. Mean field treatment 
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In Fig. 2 we show the HGPE solitons for different values of 
V/J at time t = 20/ J. We compare the continuum solu- 
tion (black line) to the solution obtained on the lattice (dashed 
lines). As can be seen, for V/J = 0.95, the lattice approxi- 
mation and the continuum solution show excellent agreement, 
up to small deviations at the boundaries. For V/J = 0.75, 
however, significant deviations occur. We further analyze this 
behavior in Fig. 3. The top panel shows the difference of the 
lattice solution to the continuum solution at t = 20/ J for a 
system of L = 40 sites as a function of V/ J. As can be seen, 
the difference is significant for all values of V/ J ^ 0.8. Only 
at such large values of V/ J the agreement is within a few per- 
cent. 

This discrepancy between the lattice and the continuum so- 
lution is to be expected: the continuum model is an approx- 
imation to the lattice model and its validity will break down 
when the size of features (e.g., the width of the soliton) of the 
analytic continuum solutions becomes comparable to the lat- 
tice spacings. This breakdown of validity can be understood in 
terms of the emission of Bogoliubov quasi-particles [55] for 
solitons which are too narrow: analogous to the excitations 
in a dilute bose gas, the Bogoliubov dispersion spectrum [51] 
shows that a narrow perturbation excites high energy modes. 
We further analyze this in Figs. 3 (c) and (d). Quasi-particles 
are emitted in the course of the time evolution, and due to 
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In this section, we compare the mean-field treatment of the 
soliton dynamics on a lattice [governed by Eqs. (5)], to the 
dynamics in the continuum [Eq. (6)]. More specifically, we 
apply the equations of motion (5) to an initial state given by 
Eqs. (7) and (9) on a finite lattice, compare to the continuum 
solution and compute the difference in the local observables 



FIG. 3. (Color online) Difference of (a) the local density p and (b) 
the condensate density p^ between the continuum evolution and the 
mean-field evolution on a lattice of L = 40 sites at times t = 20/ J 
as a function of V/ J. (c) and (d): lattice mean-field evolution of the 
local density for a broad soliton (c) with V / J — 0.65 and a narrow 
soliton (d) with V/J = 0.45. 



momentum conservation the soliton gets a velocity in the op- 
posite direction so that it starts to move away from the original 
position. The narrower the soliton, the stronger the emission 
of quasi-particles, and - as expected - the lattice approxima- 
tion becomes more and more unstable as the width of the soli- 
ton decreases, i.e., with decreasing the value of Vj J. 

Note that this behavior is reminiscent of the mechanism 
which leads to the 'light-cone' effect in correlation functions 
following a quantum quench [56-60]. In this case, the quench 
creates entangled quasi-particles on each lattice site which 
then move ballistically through the system and lead to a lin- 
ear signature in the time evolution of correlation functions. In 
this way, the velocity of the quasi-particle excitations can be 
obtained [58-60]. In a similar way, we propose that the linear 
signatures in Fig. 3 can be used to further analyze the proper- 
ties of the quasi-particles. However, this lies beyond the scope 
of the present paper so that we leave this issue open for future 
investigations. 

Due to the necessity of having a width of the soliton larger 
than a few lattice spacings, we find that we need to investi- 
gate systems with L > 30 lattice sites. In the following we 
therefore compare the lattice mean-field solution to the full 
quantum dynamics obtained by the DMRG for systems with 
L = 40 and L = 100 lattice sites and V/J > 0.9. 



B. Details for the DMRG 

We apply the adaptive t-DMRG [43, 44] for systems with 
up to L = 100 lattice sites with open boundary conditions 
(OBC). The reason we solely use OBC is that the DMRG per- 
forms far better than in the case of periodic boundary condi- 
tions, so that we can treat larger system sizes. However, at this 
point it becomes necessary to discuss the effect of the bound- 
aries: we choose system sizes and initial widths of the soli- 
tons so that there is a wide region between the soliton and the 
boundary which can be considered to be 'empty'. In Fig. 4, 
we compare the initial state for a system with L = 40 and 
L = 100 sites. As can be seen, the effect of the boundaries 
on the soliton is completely negligible. This remains so on 
time scales on which perturbations either from the boundaries 
reach the soliton or from the soliton reach the boundaries. At 
these instants of time, we stop the evolution and consider this 
to be the maximal reachable time for the given system size. 
We find that already for systems as small as L = 40 sites, the 
maximal reachable time ist > 20/ J, so that we conclude that 
the analysis which we present in the following is not affected 
by boundary effects. 

We work with the 5 = 1/2 spin system [Eq. (1)] and engi- 
neer the initial state on the lattice by imprinting a phase and 
density profile by applying an external magnetic field. More 
specifically, for the initial state we compute the ground state 
of 
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FIG. 4. (Color online) Local particle density p (blue *) and conden- 
sate density (red +) obtained by DMRG (symbols) and by the 
mean-field ansatz (solid line) for a stationary bright soliton (v = 0) 
at times t = (left panels) and at times t — 20/ J (right panels). The 
plots show results for lattice sizes of L = 100 sites (top) and L = 40 
sites (bottom). The parameters are V / J — 0.95 and = 0.45. 



with h a large multiplicative factor (~ 100) and 
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For the treatment of the dynamics of the system, we apply a 
Krylov-space variant [72] of the adaptive time dependent den- 
sity matrix renormalization group method. During the evolu- 
tion, we keep up to 1000 density-matrix eigenstates for sys- 
tems with up to L = 100 sites. We apply a time step of 
At = 0.05, resulting in a discarded weig ht of < 10-^ at 
the end of the time evolution. We estimate the error bars at 
the end of the time evolution to be smaller than the size of the 
symbols. 



IV. FULL QUANTUM DYNAMICS 

A. Soliton Stability 

The initial soliton state is prepared as discussed in Sec- 
tion IIIB and is propagated with the XXZ spin- 1/2 Hamil- 
tonian Eq. (1). Snapshots of the resulting time evolution for 
the density profile of both, the bright and the dark soliton, 
with speed ^ = are shown in Fig. 5. Since our results for 
moving solitons (v > 0) are similar, we restrict in the follow- 
ing to the case of static solitons. While the mean-field solu- 
tion remains essentially unchanged in time, the full quantum 
evolution shows some deformation of the initial state: in the 
course of the evolution, the total density profile widens as the 
peak decreases. The amount of change depends on the param- 
eters V/ J and V, and is different for the bright and the dark 
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soliton. However, as further discussed below, for V j J close 
enough to unity the difference between the quantum solution 
and the initial state remains below a few percent on a time 
scale t ~ 20/ J, where the hopping amplitude due to the map- 
ping from the spin system is J/ 2. This has to be compared to 
time scales reachable by experiments on optical lattices. For 
typical lattice depths in which a tight binding description is 
valid, the tunneling rate varies from 0.1 — 1 kHz, while the 
typical time scale of the experiments is on the order of 1-100 
milliseconds. We therefore conclude that the density profile 
suggests a stable soliton on the experimentally accessible time 
scale in the full quantum evolution. Now we turn to the con- 
densate density. Here, at t = 20/ J, the deviation from the 
mean field solution is larger. Nevertheless, as shown in Fig. 5, 
the change remains within a few percent for V/t = 0.95, so 
that we conclude that both quantities identify a stable soliton 
solution on this time scale. 

To obtain a better measure for the life time of the solitons, 
we analyze in Fig. 6 for the local observables (density p and 
condensate density , respectively) the discrepancy between 
the mean field solution and the t-DMRG evolution computed 
using Eq. (10) for different parameters. As shown in Fig. 6, 
it decreases significantly as Vj J approaches unity or as the 
speed of the soliton v (in units of the speed of sound) in- 
creases. This is associated to a widening of the initial density 
profile when increasing Vj J and a reduction of the peak am- 
plitude for larger v, so that we conclude from this analysis that 
for a variety of initial conditions the GPE and HGPE solitons 
can survive quantum fluctuations on the time scales treated. 
This is further confirmed in the following by the behavior of 
the entanglement entropy and the correlation functions. 
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FIG. 5. (Color online) Local particle density p (top) and condensate 
density p^ (bottom) obtained by the DMRG (symbols) and by the 
mean-field (solid line) propagation of the stationary — 0) bright 
(red +) and dark (blue *) solitons for times t = and 20/ J for a 
system with L = 40 sites. The parameters are V j J — 0.95 and 
po = 0.25. 
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FIG. 6. (Color online) Difference between the DMRG results and the 
lattice mean-field results for the total density (left) and for the con- 
densate density (right) for a stationary (top) and for moving (bottom) 
bright solitons (po = 0.25) on a system with L = 40 sites. 

B. Entanglement entropy and nearest neighbor correlations 

A quantity that reveals the quantum nature of a state is 
the von Neumann or entanglement entropy in the system [45] 
which is defined as 

gvN ^ _Tr(p^logp^), (13) 

with pA the reduced density matrix of a subsystem A obtained 
by tracing out the degrees of freedom of the remaining part of 
the system B. From the Schmidt decomposition 

i 

it follows that 

with |(/)^) and |(/)^) the eigenstates of the reduced density ma- 
trix of subsystem AoxB, respectively, and the eigenvalues 
of the corresponding eigenstates. This quantity gives a mea- 
sure for the entanglement between two subsystems. Since the 
initial states are product states on the lattice, S^^ is exactly 
zero at the beginning of the time evolution since only one of 
the weights is finite with = 1 while the others are exactly 
zero. If /S^^(t) remains zero (or very small) in the course of 
the time evolution, we conclude that quantum fluctuations do 
not strongly influence the nature of the initial product state, 
and so the value of S'^^{t) gives an additional measure for 
the stability of the soliton solutions. Note that there are two 
variants of this analysis: in Refs. [32, 33], the entanglement 
entropy for a subsystem of one single site is measured with 
respect to the remainder of the system. However, within the 
DMRG framework it is easier to consider the time evolution of 
the entanglement entropy for all bipartitions of the system as 
it is automatically computed in the course of the sweeps. For 
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simplicity, and since it gives a similar measure for the stabil- 
ity of the soliton, we consider here the latter. In addition, the 
behavior of this quantity for ground states of finite spin chains 
is well known from conformal field theory [61], and the nu- 
merical values can be obtained easily from the DMRG. This 
allows us to compare the values of the entanglement entropy 
during the time evolution to the ones of the strongly corre- 
lated ground state of the system which serves as a reference 
for how strongly entangled the state has become during the 
time evolution. In Fig. 7 we show typical result for the entan- 
glement entropy in ground states of the spin system Eq. (1) 
with L = 40 sites, Vj J = 0.95 and 5'*°*^^ corresponding to 
po = 0.1, 0.25 and 0.45, respectively. As can be seen, the 
numerical value in the center of the system increases with po 
and reaches 5'^^' ^ 1.35 for po = 0.45. 

A second estimate for the strength of the entanglement 
growth is to compare to the maximal possible entanglement 
entropy in a generic spin- 1/2 chain with L sites. Consider 
a bipartition of the chain into M and L — M spins with 
M < L — M. Since the dimension of the Hilbert space of 
a chain of M spins is 2^, a maximally entangled state is ob- 
tained when all Xj = 1/2^ . This state has hence an entropy 

gvN, max ^_Y^ = M log 2. 

j 

For a system of L = 40 sites and a bipartition M = L/2 we 
therefore obtain S'^^' ^ 13.86, i.e. it is a factor of - 10 
larger than the one in the ground state for the same bipartition. 

We now compare this values to the ones reached in the 
time evolution of the solitons. In Fig. 8 we display the en- 
tanglement growth of both the dark and the bright soliton at 
Po = 0.1, 0.25 and 0.45, respectively. We obtain that the 
entanglement growth is strongest at low fillings (po = 0.1), 
and it is larger for the bright soliton than for the dark one. 
For Po = 0.45, the maximum value for the bright soliton is 
5'^^ ^ 0.4, and for the dark soliton 6'^^ ^ 0.25 - both val- 
ues are significantly smaller than the one in the corresponding 
ground state, and much smaller than the one of the maximally 
entangled state. This shows that on the time scale treated, 
the state is significantly closer to a product state than to a 
strongly correlated ground state of the same system, or than 
to a maximally entangled state. Since the entanglement is not 
negligible, quantum fluctuations play an important role for the 
characterization of the state towards the end of the considered 
time evolution, but they are not strong enough to fully destroy 
the product nature of the initial state. 

Note that the entanglement growth for the bright soliton for 
po = 0.1 is significantly larger than for po = 0.45. This is 
connected to the fact that also for the local observables the 
corresponding initial state decays much faster. The entangle- 
ment entropy can be used as a measure to compare the stabil- 
ity of the initial states at po = 0.1 and po = 0.45: it appears 
that the bright soliton at po = 0.1 is about half as stable as 
the one at po = 0.45. This is reflected in the numerical val- 
ues of Sp''^^ {t) which also show approximately a factor of two 
between the two cases. 

While the entanglement entropy at the center of the system 
shows a peak for the bright soliton, it possesses a minimum for 




5 10 15 20 25 30 35 40 



position of bipartition 1 

FIG. 7. (Color online) Entanglement entropy in ground states of the 
spin system Eq. (1) for Vj J — 0.95, L = 40 sites for values of 
*^totai corresponding to the background density po = 0.1, 0.25, and 
0.45, respectively. 

the dark soliton. This can be understood by the fact that the 
dark soliton has fewer particles at the center of the system and 
so quantum fluctuations are less pronounced. Due to particle- 
hole symmetry, the dark and bright soliton evolutions for po > 
1/2 possess the same behavior. 

The behavior of the entanglement entropy can be com- 
pared to the local spin fluctuations and correlations in the 
system. In the mean field approach, at all times the coher- 
ent spin state enforces that p^^^^^ = pi(l — pi). This re- 
lation can be expressed in terms of spin observables, leading 
to (S'f )2 + (5^)2 + (5'f)2 = 1/4 on each site, realizing a 
constraint on the local spin fluctuations. In the full quantum 
dynamics, this constraint is broken, so that the initial coherent 
state becomes modified, and entanglement is induced in the 
system [62]. 

The entanglement entropy is related to the long distance 
correlations and has been extensively studied in spin systems 
[63, 68-71]. It is therefore interesting to consider the growth 
of correlations in our system in the course of the time evolu- 
tion. For simplicity, and since they are the most relevant ones 
for experiments, we consider nearest neighbor spin correlac- 
tions {Si • — {Si) • The results shown in Fig. 9 

show similar behavior to the entropy dynamics. 

C. Gaussian initial states 

In this section, we test the stability of the HGPE soliton 
solutions to modifications of the initial state. Specifically, we 
compare the time evolution of the HGPE solitons to that of a 
Gaussian initial state 

which might be easier to implement in experiments [24, 25, 
73]. We analyze the dynamics for initial states with and with- 
out a phase shift of tt across the center in order to compare 
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FIG. 8. (Color online) Entanglement entropies as a function of the subsystem size of the stationary bright soliton (top) and dark solitons 
(bottom) at different times for V/J = 0.95. 



the evolution of an initial state with a similar shape and phase 
properties as the HGPE soliton to one which has only a simi- 
lar shape. As discussed in Sec. Ill B, the initial state is created 
via a Gaussian external field. 

The obtained results are shown in Fig. 10. As can be seen, 
the Gaussian state with a phase jump remains stable and ap- 
pears to be a very good approximation to the HGPE soliton. 
In contrast, without the phase jump, the initial wave packet 
quickly disperses. Note that due to the lattice the wave packet 
can disperse by creating two peaks moving in opposite di- 
rections. This is due to the deviation of the cos(/c) disper- 
sion of the lattice from the dispersion ~ /c^ of a free particle 
and comes into appearance if the number of particles is high 
enough. 

We conclude, therefore, that once the phase jump is imple- 
mented, it is not necessary in the experiments to implement 
initial states which have exactly the form of the HGPE soli- 
tons. 



V. EXPERIMENTAL REALIZATIONS 

In this section, we discuss possible realizations of the mod- 
els introduced in Sec. II. We start with the experimental im- 
plementation of the extended Bose Hubbard model, Eq. (2), 
and its fermionic variant. The nearest neighbor interaction 
term V can be possibly generated in bosonic or spin polar- 
ized fermionic systems via long-range electric [74] or mag- 
netic [75] dipolar interactions as discussed below or with a 



short-range interaction between atoms in higher bands of the 
lattice [76]. The hard-core constraint for bosons requires in- 
creasing the interactions so that there is a large energy dif- 
ference between states with a different number of bosons per 
site. This can be achieved by tuning the scattering length via a 
Feshbach resonance [77]. Note that in this type of implemen- 
tation, in which the spin 1/2 degrees of freedom correspond to 
sites with zero and one atom, the sign of the U and V interac- 
tion is determined by the scattering length and thus is the same 
for both. In our proposal, we need attractive V interaction, so 
also the Hubbard U will be attractive. However, note that also 
in this case it is possible to realize the hard-core constraint: 
even though the states with one or zero atom per site do not 
belong to the ground state manifold, when prepared, they are 
metastable since there is no way to dump the excess energy, at 
least when prepared in the lowest band [78]. This is since the 
bandwidth in a lattice is finite, which is known to lead also to 
repulsively bound pairs [79]. In our case, however, it prevents 
double occupancies, which for |[/| oo corresponds to the 
HCB limit. A possible way to proceed then is to prepare the 
ground state in the repulsive side of the Feshbach resonance 
and then quickly ramp the magnetic field to the attractive side 
in which the evolution takes place. The atoms now are still in 
the lowest band and need to be promoted to higher bands us- 
ing, e.g., similar techniques to the ones discussed in Ref. 80. 
Note that the requirement of populating higher bands can in- 
deed lead to an additional relaxation. In the fermionic system 
the decay to the lowest band can be blocked by filling the low- 
est band. The lifetime of bosons in higher bands on the other 
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FIG. 9. (Color online) Top: contour plot of the time evolution of 
the nearest neighbor spin correlations {Si • — (Si) • 

on the whole lattice for the stationary bright soliton (left) and the 
dark soliton (right) for po = 0.25 and V j J — 0.95. Bottom: time 
evolution of the entanglement entropy for the same parameters. 
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FIG. 10. (Color online) Comparison of the time evolution of a bright 
(top) and dark (bottom) HGPE soliton (#) to the evolution of an 
initial Gaussian state with (+) and without (*) phase imprinting. 



hand does require further investigation but at least recent ex- 
periments in 2D [80] reveal that it can be 10—100 times longer 
than the characteristic time scale for intersite tunneling. 

In a recent proposal, it is shown that the XXZ spin model 
Eq. (1) and the spinless fermion model Eq. (3) both can be re- 
alized in systems of polar molecules on optical lattices, even 
though with a long-range 1 /r^ decay of the interactions rather 
than nearest-neighbor interactions only. Two different paths 
allow the study of the soliton dynamics in such experiments: 
First, as discussed in detail in Refs. [38, 39], the spin model 
Eq. (1) can be directly implemented in the case of unit filling 
(i.e., one molecule per site of the optical lattice) by select- 
ing two rotational eigenstates of the molecules which emulate 
the two spin degrees of freedom of the 5* = 1/2 chain. The 
parameters of the system can then be tuned via external DC 
electrical and microwave fields. The second implementation 
is by populating the lattice with molecules which are all in the 
same rotational eigenstate, emulating a spin polarized system. 
Since the dipolar interaction decays quickly, we presume that 
the effect of the interactions beyond nearest neighbor on the 
soliton dynamics should be very small, so that both realiza- 
tions can be used to study the soliton dynamics. We leave 
a detailed study of the effect of the interaction terms beyond 
nearest-neighbor sites on the dynamics of the solitons to fur- 
ther studies. 



VI. SUMMARY 

We have analyzed the stability and lifetime of HGPE soli- 
tons on ID lattice systems driven by a XXZ-Hamiltonian 
which can model the behavior of bosonic atoms, fermionic 
polar molecules, spin systems, and spin-polarized itinerant 
fermions on optical lattices and in condensed matter systems. 
We compare the dynamics obtained in a mean field approxi- 
mation to the full quantum evolution obtained using the adap- 
tive t-DMRG and find that the solitons remain stable under 
the full quantum evolution on time scales t ~ 20/ J, where 
J/2 is the unit of the hopping. This is quantified by the en- 
tanglement entropy which remains smaller than the one in 
the ground state of the corresponding spin system and sig- 
nificantly smaller then the one of a maximally entangled state 
on this time scale. Similar to the findings of Refs. [32, 33], 
for longer times the soliton decays. However, given the time 
scales reachable by ongoing experiments with optical lattices, 
this should suffice to identify this effect in the lab. In addition, 
we find that imperfections in the creation of the initial state 
should be of minor importance, as long as the density profile 
and the phase jump are similar to the ones of the proposed 
soliton solutions. This is exemplified by a Gaussian initial 
state, which in the case of a phase jump shows good agree- 
ment with the soliton solution, while in the absence of the 
phase jump becomes completely unstable. Due to the tunabil- 
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ity of parameters either via Feshbach resonances for atomic 
systems or via electric and microwave fields in the case of po- 
lar molecules, the possibility of realizing both, bright and dark 
solitons in strongly interacting systems, adds a new paradigm 
to the existence of coherent non-linear modes in systems of 
ultracold quantum gases. 
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